function obj = locproj_var_sur(varargin)
% Estimates the local projections to be used for the pathwise equality
% test (a system of seemingly unrelated regressions) presented in Tables 3
% and 4. 
global PP H
    switch length(varargin)
        case 6
            y       =  varargin{1};
            x       =  varargin{2};
            w       =  varargin{3};
            h1      =  varargin{4};
            H       =  varargin{5};
            type    =  varargin{6};

        case 8
            y       =  varargin{1};
            x       =  varargin{2};
            w       =  varargin{3};
            h1      =  varargin{4};
            H       =  varargin{5};
            type    =  varargin{6};            
            r       =  varargin{7};
            lam     =  varargin{8};

        otherwise
            error('wrong number of input arguments')
    end    

    obj     =   struct();
    std     =   strcmp('reg',type);
    T       =   length(y);
    HR      =   H + 1 - h1;

    % Constructr the B-spline basis functions:
    if ~std
            B       =   bspline((h1:H)',h1,H+1,H+1-h1,3);
            K       =   size(B,2);
    else
            K       =   HR;
    end

    % Building up the regression representation of the local projection:
    idx     =   nan((H+1)*T,2);
    Y       =   nan((H+1)*T,1);
    Xb      =   zeros((H+1)*T,K);
    Xc      =   zeros((H+1)*T,HR,size(w,2)+1);
    w       =   [ones(T,1) w];

    for tx = 1:T

        idx_beg                             =   (tx-1)*HR + 1;
        idx_end                             =   tx*HR;
        idx(idx_beg:idx_end,1)              =   tx;
        idx(idx_beg:idx_end,2)              =   1:H;
        y_range                             =   (tx:min(tx+H-1,T)); 
        Y(idx_beg:idx_end)                  =   [y(y_range);nan(HR-length(y_range),1)];

        if std
                Xb(idx_beg:idx_end,:)           =   eye(HR)*x(tx);
        else
                Xb(idx_beg:idx_end,:)           =   B*x(tx);
        end

        for ix = 1:size(w,2)
                Xc(idx_beg:idx_end,:,ix)   =   eye(HR)*w(tx,ix);
        end

    end

    X   =   Xb;
        for ix = 1:size(w,2)
            X           =   [X Xc(:,:,ix)];
        end

    select          =   isfinite(Y);  
    idx             =   idx(select,:);
    Y               =   Y(select);
    X               =   X(select,:);
    Y               =   Y(1:(T-H-1)*H,:);
	X               =   X(1:(T-H-1)*H,:);
    % Estimation:
    IR              =   zeros(HR+1,1);

    if std
        theta               =   (X'*X )\(X'*Y);
        IR((1+h1):(1+HR))   =   theta(1:K);
    else
        P                   =   zeros(size(X,2));
        D                   =   eye(K);
        
        for kx = 1:r 
            D                   =   diff(D);
        end
        
        P(1:K,1:K)          =   D'*D;
        theta               =   (X'*X + lam*P)\(X'*Y);

        IR((1+h1):(1+HR))   =   B*theta(1:K);
    end

    % pack everything up
    obj.T       =   T;
    obj.HR      =   HR;
    obj.K       =   K;

    if ~std
        P2          =   P;
        obj.B       =   B;
        obj.P2      =   P2;
    end

    obj.idx     =   idx;
    obj.Y       =   Y;
    obj.X       =   X;
    obj.theta   =   theta;
    obj.IR      =   IR;
 end

function B = bspline(x, xl, xr, ndx, bdeg)
    dx      =   (xr - xl) / ndx;
    t       =   xl + dx * [-bdeg:ndx-1];
    T       =   (0 * x + 1) * t;
    X       =   x * (0 * t + 1);
    P       =   (X - T) / dx;
    B       =   (T <= X) & (X < (T + dx));
    r       =   [2:length(t) 1];
    for kx = 1:bdeg
        B       =   (P.*B + (kx+1-P).*B(:,r))/kx;
    end
end
